% Copyright (C) 2014-19 Benjamin Born and Johannes Pfeifer 
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This code is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% For a copy of the licencse,
% see <http://www.gnu.org/licenses/>.

%% Table 3 in the paper

% Set the variable of interest in line 9. 1: g forecast error (lower panel), 3: Default premium (upper panel)
% Set unconditional vs conditional model in line 11
%% Setup

clear; close all; clc
addpath('../Auxiliary_Files')
table_variable = 1; % 1: g forecast error, 3: risk premium

unconditional_model_dummy = 0;

impulse_name = 'PC_spread_group_specific';

load('../data_array_for_regression_stacked_by_variable_CDS_new')

lag_number=4;
max_irf_horizon = 2;
time_fixed_effects_dummy=0;
country_fixed_effects_dummy=1;
log_spread_dummy=0;
firstdiff_spread_dummy=0;

add_residual_dummy = 1; %adds residuals from previous step to regression to increase efficiency;

regressor_var_names={impulse_name}; %do not add spreads here, they are added automatically depending on indicator_dummy
dependent_var_names={'g_growth_fe','y_growth_realization','spread_end_quarter'}; %make sure spread is consistent with indicator if requested!

Indicator_position = pos.empirical_cdf_country_group_specific_end_quarter;
[regressor_var_names] = check_dependent_and_regressors_for_consistency(dependent_var_names,regressor_var_names,'spread_end_quarter',log_spread_dummy,firstdiff_spread_dummy);
regressor_var_names = regressor_var_names(end-1);

% No US/Germany
[row,col] = find(data_array_for_regression_stacked_by_variable(:,:,pos.country_ident)==39 ...
    | data_array_for_regression_stacked_by_variable(:,:,pos.country_ident)==12);
for ii=1:length(row)
    data_array_for_regression_stacked_by_variable(row(ii),col(ii),3:end)=NaN;
end

%% Shock Variable

fe_shocks_temp = data_array_for_regression_stacked_by_variable(:,:,pos.(impulse_name));

fe_shocks_negative = zeros(size(fe_shocks_temp));
fe_shocks_negative(fe_shocks_temp<0 | isnan(fe_shocks_temp)) = fe_shocks_temp(fe_shocks_temp<0 | isnan(fe_shocks_temp));
fe_shocks_positive = zeros(size(fe_shocks_temp));
fe_shocks_positive(fe_shocks_temp>0 | isnan(fe_shocks_temp)) = fe_shocks_temp(fe_shocks_temp>0 | isnan(fe_shocks_temp));
if any(any(fe_shocks_negative>0)) || any(any(fe_shocks_positive<0 ))
    error('Sign is wrong')
end

if unconditional_model_dummy == 1
    data_array_for_regression_impulse_only=cat(3,fe_shocks_temp);
else
    % neg_shock*F(z)
    data_array_for_regression_g_stress_neg=fe_shocks_negative.*...
        data_array_for_regression_stacked_by_variable(:,:,Indicator_position);
    % neg_shock*(1-F(z))
    data_array_for_regression_g_nostress_neg=fe_shocks_negative.*...
        (1-data_array_for_regression_stacked_by_variable(:,:,Indicator_position));
    
    % pos_shock*F(z)
    data_array_for_regression_g_stress_pos=fe_shocks_positive.*...
        data_array_for_regression_stacked_by_variable(:,:,Indicator_position);
    % pos_shock*(1-F(z))
    data_array_for_regression_g_nostress_pos=fe_shocks_positive.*...
        (1-data_array_for_regression_stacked_by_variable(:,:,Indicator_position));
    
    % shock*F(z)
    data_array_for_regression_g_stress=fe_shocks_temp.*...
        data_array_for_regression_stacked_by_variable(:,:,Indicator_position);
    % shock*(1-F(z))
    data_array_for_regression_g_nostress=fe_shocks_temp.*...
        (1-data_array_for_regression_stacked_by_variable(:,:,Indicator_position));
    
    data_array_for_regression_impulse_only=cat(3,...
        data_array_for_regression_g_stress,data_array_for_regression_g_nostress...
        );
end
data_array_for_regression = data_array_for_regression_impulse_only;

%% Lagged regressors

% ATTENTION: Spreads always last regressor

for var_iter=1:length(regressor_var_names)
    if unconditional_model_dummy
        data_array_control=data_array_for_regression_stacked_by_variable(:,:,pos.([regressor_var_names{var_iter},'_lag_1']):pos.([regressor_var_names{var_iter},'_lag_',num2str(lag_number)]));
        data_array_for_regression=cat(3,data_array_for_regression,data_array_control);
    else
        %F(z)*X_{t-k}
        data_array_state_1 = data_array_for_regression_stacked_by_variable(:,:,pos.([regressor_var_names{var_iter},'_lag_1']):pos.([regressor_var_names{var_iter},'_lag_',num2str(lag_number)])).*...
            (repmat(data_array_for_regression_stacked_by_variable(:,:,Indicator_position),[1,1,lag_number]));
        %(1-F(z))*X_{t-k}
        data_array_state_2 = data_array_for_regression_stacked_by_variable(:,:,pos.([regressor_var_names{var_iter},'_lag_1']):pos.([regressor_var_names{var_iter},'_lag_',num2str(lag_number)])).*...
            (repmat(1-data_array_for_regression_stacked_by_variable(:,:,Indicator_position),[1,1,lag_number]));
        
        data_array_for_regression=cat(3,data_array_for_regression,data_array_state_1,data_array_state_2);
    end
end

%% Run regression

if unconditional_model_dummy == 1
    theta_mat=NaN(1+max_irf_horizon,1,length(dependent_var_names));
    se_mat=NaN(1+max_irf_horizon,1,length(dependent_var_names));
else
    theta_mat=NaN(1+max_irf_horizon,2,length(dependent_var_names));
    se_mat = NaN(1+max_irf_horizon,2,length(dependent_var_names));
    test_stat_mat = NaN(1+max_irf_horizon,length(dependent_var_names));
    test_p_mat = NaN(1+max_irf_horizon,length(dependent_var_names));
end

for var_iter=1:length(dependent_var_names)
    for horizon_iter = 0:max_irf_horizon
        %construct matrices without NaN
        if horizon_iter==0
            dependent_variable=squeeze(data_array_for_regression_stacked_by_variable(:,:,pos.(dependent_var_names{var_iter})));
            yhat_full=NaN(size(dependent_variable));
            [dependent_variable_for_run,data_array_for_regression_for_run,time_indices_non_NaN,country_indices_non_NaN]=create_regression_matrices_no_NaN(dependent_variable,data_array_for_regression,data_array_for_regression_stacked_by_variable,pos,country_indicator_names_mapping,time_fixed_effects_dummy,country_fixed_effects_dummy);
            %if var_iter==1 && horizon_iter==0
            %    generate_descriptives_table;
            %end
            if var_iter==1 && horizon_iter==0
                generate_nobs_ncountry;
            end
            %run regression
            [theta,stdDK,~,CovDK,yhat] = HszDk5cPs(dependent_variable_for_run,ones(size(dependent_variable_for_run,1),1),data_array_for_regression_for_run,1,8,1);
            
            if add_residual_dummy && ~strcmp(impulse_name,dependent_var_names(var_iter)) %when regression G on itself, residuals are 0
                yhat_full(time_indices_non_NaN,country_indices_non_NaN)=yhat;
                resids=dependent_variable-yhat_full;
            else
                resids=[];
            end
        else
            dependent_variable=squeeze(data_array_for_regression_stacked_by_variable(:,:,pos.([dependent_var_names{var_iter},'_lead_',num2str(horizon_iter)])));
            yhat_full=NaN(size(dependent_variable));
            
            [dependent_variable_for_run,data_array_for_regression_for_run,time_indices_non_NaN,country_indices_non_NaN]=create_regression_matrices_no_NaN(dependent_variable,cat(3,data_array_for_regression,resids),data_array_for_regression_stacked_by_variable,pos,country_indicator_names_mapping,time_fixed_effects_dummy,country_fixed_effects_dummy);
            %run regression
            [theta,stdDK,~,CovDK,yhat] = HszDk5cPs(dependent_variable_for_run,ones(size(dependent_variable_for_run,1),1),data_array_for_regression_for_run,1,8,1);
            if add_residual_dummy
                yhat_full(time_indices_non_NaN,country_indices_non_NaN)=yhat;
                resids=dependent_variable-yhat_full;
            else
                resids=[];
            end
        end
        
        %save coefficients
        if unconditional_model_dummy
            theta_mat(1+horizon_iter,:,var_iter) = theta(1,1);
            se_mat(1+horizon_iter,:,var_iter)  = stdDK(1,1)';
        else unconditional_model_dummy == 0
            theta_mat(1+horizon_iter,1:2,var_iter) = theta(1:2,1)';
            se_mat(1+horizon_iter,:,var_iter)  = stdDK(1:2,1)';
            Rmat = zeros(1,length(theta));
            Rmat(1,1:2) = [1,-1];
            Wstat = (Rmat*theta)'/(Rmat*CovDK*Rmat')*(Rmat*theta);
            test_stat_mat(1+horizon_iter, var_iter) = Wstat;
            test_p_mat(1+horizon_iter, var_iter) = 1 - chi2cdf(Wstat,1);
        end
    end
end

%% Output table

scaling = 100;

% only keep contemporaneous spread part

theta_mat_spread = squeeze(theta_mat(1:2,:,table_variable));
se_mat_spread    = squeeze(se_mat(1:2,:,table_variable));
p_mat_spread     = (1-normcdf(abs(theta_mat_spread)./se_mat_spread))*2;

% for syntax, see latexTable example in Data/Data/Auxiliary_Files
input.transposeTable = 0;
input.dataNanString = '-';
input.tableColumnAlignment = 'l';
input.tableBorders = 0;
input.booktabs = 1;

if unconditional_model_dummy
    input.data = [scaling .* theta_mat_spread; abs(scaling) .* se_mat_spread; p_mat_spread; n_countries; n_obs];
    input.tableRowLabels = {'coeff. 0','coeff. 1','se 0','se 1','pValue 0','pValue 1','n_country' 'n_obs'};
    input.tableColLabels = {'Linear'};
    input.dataFormat = {'%2.3f',1};
    input.tableCaption = 'Unconditional';
    input.tableLabel = 'uncond';
else unconditional_model_dummy
    input.data = [scaling .* theta_mat_spread; abs(scaling) .* se_mat_spread; p_mat_spread; [n_countries n_obs]];
    input.tableRowLabels = {'coeff. 0','coeff. 1','se 0','se 1','pValue 0','pValue 1','n_obs_country'};
    input.tableColLabels = {'Stress','NoStress'};
    input.dataFormat = {'%2.3f',2};
    input.tableCaption = 'Conditional';
    input.tableLabel = 'cond';
end

input.makeCompleteLatexDocument = 0;
latex = latexTable(input);
